options(
java.parameters = c("-XX:+UseConcMarkSweepGC", "-Xmx36g")
)
gc()Análisis del paisaje
El paquete usado es {rflsgen}, Neutral Landscape Generator with Targets on Landscape Indices.
Antes que nada hay que ejecutar los siguientes comandos:
Esto permite aumentar la memoria disponible para las funciones de {rflsgen}.
Los paquetes usados son los siguientes:
library(tidyverse)
library(gt)
library(landscapemetrics)
library(NLMR)
library(rflsgen)
library(terra)1 Lectura de datos
Inicialmente hay que leer los ráster de uso de suelo (r) y modelo digital de elevación (dem).
Código
r <- rast("datos/r.tif")
dem <- rast("datos/dem.tif")Las clases de r son:
| Valor | Tipo |
|---|---|
| 0 | No data |
| 10 | Forest |
| 20 | Scrubland |
| 30 | Grassland |
| 40 | Agriculture |
| 50 | Urban |
| 60 | Sparse vegetation - bare soil |
| 80 | Water body |
Agrego las categorías a r e incorporo colores. Visualizo ambos ráster.
Código
d <- read_delim(
file = "datos/LULCcode.txt",
delim = "$ ",
col_names = c("value"),
show_col_types = FALSE
) |>
separate_wider_delim(
delim = " ",
cols = value,
names = c("value", "tipo"),
too_many = "merge"
) |>
mutate(value = as.integer(value)) |>
filter(value != 0) |>
mutate(col = c(
"#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#A65628", "#F781BF"
))
col_df <- d |>
select(-tipo) |>
as.data.frame()
# incorporo categorías y colores
levels(r) <- d
coltab(r) <- col_df2 Extracción de la estructura del paisaje
Se seleccionan las clases de interés del ráster de cobertura del suelo.
Código
clases <- c(30, 40, 50) # Grassland, Agriculture, Urban
estructura <- flsgen_extract_structure_from_raster(
raster_file = r,
focal_classes = clases,
connectivity = 8
)A partir de estructura y dem se genera un paisaje, indicando dependencia del terreno, conectividad, rugosidad y parámetros de posición.
Para visualizar el efecto de la dependencia del terreno (terrain_dependency) se generan 4 escenarios con valores crecientes: 0.1, 0.3, 0.7 y 1.
Código
f_conectividad <- function(x) {
print(paste0("Conectividad: ", x))
r_land <- flsgen_generate(
estructura,
terrain_file = dem,
terrain_dependency = x,
connectivity = 8,
epsg = "EPSG:22174",
resolution_x = res(r)[1],
resolution_y = res(r)[2],
x = ext(dem)$"xmin",
y = ext(dem)$"ymax"
)
writeRaster(r_land, filename = paste0("rasters/land-", x, ".tif"))
}
conectividad <- c(.1, .3, .7, 1)
walk(conectividad, f_conectividad)Los escenarios generados son:
El paquete {rflsgen} genera múltiples escenarios pero manteniendo constantes las métricas del paisaje.
Para el cálculo de las métricas se utiliza el paquete {landscapemetrics}. Se calcularon: Mean of patch area y Number of patches.
Código
f_area_mn <- function(x) {
lsm_c_area_mn(land_raster[[x]], directions = 8) |>
mutate(conectividad = conectividad[x]) |>
mutate(metrica = "Mean of patch area")
}
df_area_mn <- map(1:length(conectividad), f_area_mn) |>
list_rbind()
write_tsv(df_area_mn, "indices/df_area_mn.tsv")
f_np <- function(x) {
lsm_c_np(land_raster[[x]], directions = 8) |>
mutate(conectividad = conectividad[x]) |>
mutate(metrica = "Number of patches")
}
df_np <- map(1:length(conectividad), f_np) |>
list_rbind()
write_tsv(df_np, "indices/df_np.tsv")Los valores se obtienen para cada clase elegida previamente. Se observa que para cada valor de conectividad las métricas coinciden.
Código
read_tsv("indices/df_area_mn.tsv", show_col_types = FALSE) |>
filter(class != -1) |>
mutate(value = round(value, 2)) |>
inner_join(d_r, by = join_by(class == value)) |>
pivot_wider(
names_from = conectividad,
values_from = value,
names_prefix = "Dependencia: "
) |>
select(Tipo = tipo, starts_with("Dependencia")) |>
gt() |>
tab_header("Mean of patch area") |>
tab_style(
locations = cells_body(columns = starts_with("Dependencia")),
style = cell_text(font = "JetBrains Mono", color = "black")
) |>
fmt_number(
columns = everything(),
dec_mark = ",",
decimals = 2,
sep_mark = "."
)| Mean of patch area | ||||
|---|---|---|---|---|
| Tipo | Dependencia: 0.1 | Dependencia: 0.3 | Dependencia: 0.7 | Dependencia: 1 |
| Grassland | 3,52 | 3,52 | 3,52 | 3,52 |
| Agriculture | 29,54 | 29,54 | 29,54 | 29,54 |
| Urban | 1,88 | 1,88 | 1,88 | 1,88 |
Código
read_tsv("indices/df_np.tsv", show_col_types = FALSE) |>
filter(class != -1) |>
mutate(value = round(value, 2)) |>
inner_join(d_r, by = join_by(class == value)) |>
pivot_wider(
names_from = conectividad,
values_from = value,
names_prefix = "Dependencia: "
) |>
select(Tipo = tipo, starts_with("Dependencia")) |>
gt() |>
tab_header("Number of patches") |>
tab_style(
locations = cells_body(columns = starts_with("Dependencia")),
style = cell_text(font = "JetBrains Mono", color = "black")
) |>
fmt_number(
columns = everything(),
dec_mark = ",",
decimals = 0,
sep_mark = "."
)| Number of patches | ||||
|---|---|---|---|---|
| Tipo | Dependencia: 0.1 | Dependencia: 0.3 | Dependencia: 0.7 | Dependencia: 1 |
| Grassland | 45.871 | 45.871 | 45.871 | 45.871 |
| Agriculture | 12.308 | 12.308 | 12.308 | 12.308 |
| Urban | 15.311 | 15.311 | 15.311 | 15.311 |
En la generación de los escenarios pueden suceder dos errores:
Falta de memoria, al no contar con la suficiente capacidad de procesamiento. Para solucionarlo se puede reducir el tamaño de la región de interés, aumentar el tamaño de píxel o disminuir el número de clases.
Incapacidad de generar el ráster, mostrando el siguiente mensaje:
Could not generate a raster satisfying the input landscape structureEsto puede solucionarse eligiendo únicamente las clases de interés.
3 Estructura del paisaje propia
Genero el caso I y el caso II, cada uno con dos clases y métricas del paisaje diferentes.
Luego, para verificar los escenarios generados, calculo las métricas y verifico que se cumplan las condiciones dadas.
Defino dos clases y sus métricas.
Código
tibble(
Clase = c("Clase 1", "Clase 2"),
`Number of patches` = c("[700, 800]", "[800, 1100]"),
`Patch area` = c("[5500, 14600]", "[9500, 11060]"),
`Proportion of landscape` = c("10%", "-")
) |>
gt() |>
cols_label(
`Number of patches` = md("Number of<br>patches"),
`Proportion of landscape` = md("Proportion of<br>landscape"),
"Clase" = ""
) |>
tab_style(
locations = cells_body(columns = -Clase),
style = cell_text(font = "JetBrains Mono", color = "black")
) |>
tab_header(title = md("**Caso I**"))| Caso I | |||
|---|---|---|---|
| Number of patches |
Patch area | Proportion of landscape |
|
| Clase 1 | [700, 800] | [5500, 14600] | 10% |
| Clase 2 | [800, 1100] | [9500, 11060] | - |
Código
tibble(
Clase = c("Clase A", "Clase B"),
`Number of patches` = c("[700, 800]", "[700, 800]"),
`Patch area` = c("[15500, 114600]", "[15500, 114600]"),
`Proportion of landscape` = c("-", "25%")
) |>
gt() |>
cols_label(
`Number of patches` = md("Number of<br>patches"),
`Proportion of landscape` = md("Proportion of<br>landscape"),
"Clase" = ""
) |>
tab_style(
locations = cells_body(columns = -Clase),
style = cell_text(font = "JetBrains Mono", color = "black")
) |>
tab_header(title = md("**Caso II**"))| Caso II | |||
|---|---|---|---|
| Number of patches |
Patch area | Proportion of landscape |
|
| Clase A | [700, 800] | [15500, 114600] | - |
| Clase B | [700, 800] | [15500, 114600] | 25% |
# CASO I
cls_1 <- flsgen_create_class_targets(
class_name = "Clase 1",
NP = c(700, 800),
AREA = c(5500, 14600),
PLAND = c(10, 10)
)
cls_2 <- flsgen_create_class_targets(
class_name = "Clase 2",
NP = c(800, 1100),
AREA = c(9500, 11060)
)
# CASO II
cls_A <- flsgen_create_class_targets(
class_name = "Clase A",
NP = c(700, 800),
AREA = c(15500, 114600)
)
cls_B <- flsgen_create_class_targets(
class_name = "Clase B",
NP = c(700, 800),
AREA = c(15500, 114600)
)Las métricas y los valores tienen que ser coherentes y no pueden adoptar cualquier número. En caso de indicar valores incorrectos, se muestra el siguiente mensaje de error:
User targets cannot be satisfiedCreo un conjunto con las clases y obtengo la estructura, utilizando el dem como máscara.
Código
# CASO I
objetivos <- flsgen_create_landscape_targets(
mask_raster = dem,
classes = list(cls_1, cls_2)
)
estructura_custom <- flsgen_structure(objetivos)
# CASO II
objetivos2 <- flsgen_create_landscape_targets(
mask_raster = dem,
classes = list(cls_A, cls_B)
)
estructura_custom2 <- flsgen_structure(objetivos2)Genero un ráster a partir de las clases creadas y su estructura ubicados sobre el dem original.
Código
r_custom <- flsgen_generate(
estructura_custom,
terrain_file = dem,
terrain_dependency = 1,
connectivity = 8,
epsg = "EPSG:22174",
resolution_x = res(r)[1],
resolution_y = res(r)[2],
x = ext(dem)$"xmin",
y = ext(dem)$"ymax"
)
writeRaster(r_custom, "rasters/r_custom.tif", overwrite = TRUE)
r_custom2 <- flsgen_generate(
estructura_custom2,
terrain_file = dem,
terrain_dependency = 1,
connectivity = 8,
epsg = "EPSG:22174",
resolution_x = res(r)[1],
resolution_y = res(r)[2],
x = ext(dem)$"xmin",
y = ext(dem)$"ymax"
)
writeRaster(r_custom2, "rasters/r_custom2.tif", overwrite = TRUE)Visualizo los resultados de ambos casos.
Código
r_custom <- rast("rasters/r_custom.tif")
d_custom <- tibble(
value = as.integer(c(0, 1)),
tipo = c("Clase 1", "Clase 2"),
col = c("#E41A1C", "#4DAF4A")
) |>
add_row(
value = as.integer(-1),
tipo = "X",
col = "grey95"
)
col_df_custom <- d_custom |>
select(-tipo) |>
as.data.frame()
# incorporo categorías y colores
levels(r_custom) <- d_custom
coltab(r_custom) <- col_df_custom
plot(r_custom, axes = FALSE, mar = c(1, 1, 1, 4), main = "Caso I")
r_custom2 <- rast("rasters/r_custom2.tif")
d_custom2 <- tibble(
value = as.integer(c(0, 1)),
tipo = c("Clase A", "Clase B"),
col = c("#E41A1C", "#4DAF4A")
) |>
add_row(
value = as.integer(-1),
tipo = "X",
col = "grey95"
)
col_df_custom2 <- d_custom2 |>
select(-tipo) |>
as.data.frame()
# incorporo categorías y colores
levels(r_custom2) <- d_custom2
coltab(r_custom2) <- col_df_custom2
plot(r_custom2, axes = FALSE, mar = c(1, 1, 1, 4), main = "Caso II")Verifico que el ráster creado tenga las métricas previamente definidas para cada clase.
Código
read_tsv("datos/df_custom.tsv", show_col_types = FALSE) |>
filter(class != -1) |>
mutate(
clase = if_else(class == 0, "Clase 1", "Clase 2")
) |>
pivot_wider(
names_from = metrica,
values_from = value,
id_cols = clase
) |>
select(clase, `Number of patches`, `Percentage of landscape of class`) |>
gt() |>
cols_label(
`Number of patches` = md("Number of<br>patches"),
`Percentage of landscape of class` = md("Percentage of<br>landscape of class"),
"clase" = ""
) |>
tab_style(
locations = cells_body(columns = -clase),
style = cell_text(font = "JetBrains Mono", color = "black")
) |>
fmt_number(
columns = 2,
dec_mark = ",",
decimals = 0,
sep_mark = "."
) |>
fmt_number(
columns = 3,
dec_mark = ",",
decimals = 1,
sep_mark = "."
) |>
tab_header(title = md("**Caso I**"))| Caso I | ||
|---|---|---|
| Number of patches |
Percentage of landscape of class |
|
| Clase 1 | 791 | 10,0 |
| Clase 2 | 806 | 11,5 |
Código
read_tsv("datos/df_custom2.tsv", show_col_types = FALSE) |>
filter(class != -1) |>
mutate(
clase = if_else(class == 0, "Clase A", "Clase B")
) |>
pivot_wider(
names_from = metrica,
values_from = value,
id_cols = clase
) |>
select(clase, `Number of patches`, `Percentage of landscape of class`) |>
gt() |>
cols_label(
`Number of patches` = md("Number of<br>patches"),
`Percentage of landscape of class` = md("Percentage of<br>landscape of class"),
"clase" = ""
) |>
tab_style(
locations = cells_body(columns = -clase),
style = cell_text(font = "JetBrains Mono", color = "black")
) |>
fmt_number(
columns = 2,
dec_mark = ",",
decimals = 0,
sep_mark = "."
) |>
fmt_number(
columns = 3,
dec_mark = ",",
decimals = 1,
sep_mark = "."
) |>
tab_header(title = md("**Caso II**"))| Caso II | ||
|---|---|---|
| Number of patches |
Percentage of landscape of class |
|
| Clase A | 701 | 16,3 |
| Clase B | 753 | 25,0 |